!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!! This module summarizes the attributes of the master element (also
!! said reference element) for a simplex based \f$d\f$-dimensional
!! grid.
!!
!! \n
!!
!! The master element is the \f$d\f$-simplex \f$\hat{E}\f$ of
!! \f$\mathbb{R}^d\f$ with vertices \f$v^E_1=(0,0,\ldots,0)\f$,
!! \f$v^E_2=(1,0,\ldots,0), \ldots, v^E_{d+1}=(0,0,\ldots,1)\f$,
!! in this order (which also implies positive orientation). The sides
!! of \f$\hat{E}\f$ are denoted by \f$\hat{F}\f$, and for each side an
!! ordering of the \f$d\f$ vertices is specified such that the side
!! normal \f$n^F\f$ is directed outwards. For more details concerning
!! the structure of a \f$d\f$-dimensional grid see \c mod_grid.
!<----------------------------------------------------------------------
module mod_master_el

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_octave_io, only: &
   mod_octave_io_initialized, &
   write_octave

 use mod_octave_io_perms, only: &
   mod_octave_io_perms_initialized, &
   write_octave

 use mod_perms, only: &
   mod_perms_initialized, &
   t_perm,      &
   operator(*), &
   perm_table,  &
   idx,         &
   fact,        &
   comb_table

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_master_el_constructor, &
   mod_master_el_destructor,  &
   mod_master_el_initialized, &
   t_me, new_me, clear,   &
   t_lagnodes, lag_nodes, &
   add_nodes,             &
   write_octave

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members

 !> a list of simplexes, represented as a two-dimensional integer
 !! array: <tt>s(:,i)</tt> are the vertices of the <tt>i</tt>-th
 !! simplex.
 !! \bug the two fields \c d and \c ns should be type parameters.
 !<
 type t_simpl_list
   integer :: d  !< dimension
   integer :: ns !< number of simplexes
   integer, allocatable :: s(:,:) !< vertices (d+1,ns)
 end type t_simpl_list

 !> master element
 !> \bug The field \c d should be a derived type parameter
 type t_me
   integer :: d !< dimension
   !> <tt>isv(:,is)</tt> contains the nodes of the side \c is of
   !! master element for \c is between 1 and <tt>d+1</tt>.
   !<
   integer, allocatable :: isv(:,:)
   !> <tt>children(i)</tt> is a t_simpl_list object containing all the
   !! boundary \f$i\f$-simplexes of the master element, for
   !! <tt>i=0,d-1</tt>.
   !! \note The simplexes are stored in increasing lexicographic
   !! order, since this is the only ordering which makes sense for any
   !! simplex dimension. This means that the sides, which appear as
   !! \f$(d-1)\f$-simplexes in \c children as well as in \c isv, have
   !! two different orderings and orientations in the two fields. In
   !! \c children, the sides are not ordered according to the opposite
   !! vertex and their normal is not necessarily pointing outwards.
   !<
   type(t_simpl_list), allocatable :: children(:)
   !> \c pi_tab lists the \f$d!\f$ permutations of order \f$d\f$ in
   !! lexicographic order.
   !<
   type(t_perm), allocatable :: pi_tab(:)
   real(wp), allocatable :: xv(:,:) !< vertex coordinates
   !> <tt>n(:,is)</tt> is the outward normal of the <tt>is</tt>-th
   !! side.
   !<
   real(wp), allocatable :: n(:,:)
   real(wp), allocatable :: a(:) !< side areas
   real(wp) :: vol !< volume
   real(wp) :: voldm1 !< volume of the reference \f$(d-1)\f$-simplex
 end type t_me

 !> Lagrangian nodes.
 !!
 !! This type collects the Lagrangian nodes for a \f$d\f$-dimensional
 !! simplex. The field \c stab is a permutation table giving the node
 !! permutation corresponding to each vertex permutation (this
 !! information is analogous to the field \c stab of a \c t_base
 !! object). 
 !!
 !! In a typical use, only internal nodes are stored in a \c
 !! t_lagnodes object, because it is assumed that boundary nodes will
 !! be stored in another lower dimensional \c t_lagnodes object.
 !! Hence, the entire set of Lagrangian nodes for a \f$d\f$-simplex
 !! will be stored in an array \code
 !!   type(t_lagnodes) :: lagnodes(1:d)
 !! \endcode
 !! with <tt>lagnodes(i)</tt> giving the nodes on \f$i\f$-dimensional
 !! sub-simplexes of the \f$d\f$-simplex. \f$0\f$-simplexes coincide
 !! with vertices, and they should be handled separately.
 !! \bug The fields <tt>d k nn</tt> should be derived type parameters
 !<
 type t_lagnodes
   integer :: d !< dimension
   integer :: nn !< number of nodes
   real(wp), allocatable :: x(:,:) !< nodes (d,nn)
   integer, allocatable :: stab(:,:) !< permutation table ((d+1)!,nn)
 end type t_lagnodes

! Module variables

 ! public members
 logical, protected ::               &
   mod_master_el_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_master_el'

 interface clear
   module procedure clear_me
 end interface

 interface write_octave
   module procedure write_me, write_simpl_list, write_lagnodes
 end interface write_octave

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_master_el_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or.  &
       (mod_kinds_initialized.eqv..false.) .or.     &
       (mod_octave_io_initialized.eqv..false.) .or. &
       (mod_octave_io_perms_initialized.eqv..false.) .or. &
       (mod_perms_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_master_el_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_master_el_initialized = .true.
 end subroutine mod_master_el_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_master_el_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_master_el_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_master_el_initialized = .false.
 end subroutine mod_master_el_destructor

!-----------------------------------------------------------------------
 
 !> Define a t_me object
 pure subroutine new_me(me,d)
  integer, intent(in) :: d
  type(t_me), intent(out) :: me
  character(len=*), parameter :: &
    this_sub_name = 'new_me'

  integer :: is, i

   me%d = d
   allocate( me%isv(d,d+1) )
   do is=1,d+1 ! side loop
     me%isv(:,is) = mod( is+(/(i, i=1,d)/)-1 , d+1 )+1
     ! the following parity change works empirically up to d=20;
     ! formal proof still lacking...
     if((-1)**(d*is).eq.-1) then
       i = me%isv(1,is)
       me%isv(1,is) = me%isv(d,is)
       me%isv(d,is) = i
     endif
   enddo
   allocate( me%children(0:d-1) )
   do i=0,d-1
     call comb_table(me%children(i)%s,(/(is, is=1,d+1)/),i+1)
     me%children(i)%d = i
     me%children(i)%ns = size(me%children(i)%s,2)
   enddo
   call perm_table(me%pi_tab,d) ! side vertex permutations
 
   ! geometry
   allocate( me%xv(d,d+1) , me%n(d,d+1) , me%a(d+1) )
   me%xv = 0.0_wp
   do is=2,d+1
     me%xv(is-1,is) = 1.0_wp
   enddo
   me%n = 0.0_wp
   me%n(:,1) = 1.0_wp/sqrt(real(d,wp))
   do is=2,d+1
     me%n(is-1,is) = -1.0_wp
   enddo
   me%a(1) = sqrt(real(d,wp))/real(fact(d-1),wp)
   me%a(2:d+1) = 1.0_wp/real(fact(d-1),wp)
   me%vol      = 1.0_wp/real(fact(d)  ,wp)
   me%voldm1   = 1.0_wp/real(fact(d-1),wp)

 end subroutine new_me
 
!-----------------------------------------------------------------------
 
 !> Deallocate a master element object
 pure subroutine clear_me(me)
  type(t_me), intent(out) :: me
 end subroutine clear_me

!-----------------------------------------------------------------------

 !> Compute the Lagrangian nodes for a basis of the space
 !! \f$\mathbb{P}_k(\mathbb{R}^d)\f$ on \f$\hat{K}\f$. 
 !!
 !! For a \f$d\f$-simplex, the Lagrangian nodes can be associated with
 !! any of the \f$e\f$-simplexes representing the \f$e\f$-faces, for
 !! \f$e=0,\ldots,d-1\f$, as well as with the \f$d\f$-simplex itself.
 !! This means that there are nodes coinciding with the vertices,
 !! nodes placed on the 1-dimensional edges and so on up to nodes
 !! placed on the \f$(d-1)\f$-dimensional sides and finally internal
 !! nodes. To obtain a conforming finite element space, all the nodes
 !! placed on the boundary must be symmetric (i.e. the set of the
 !! nodes must be closed with respect to all the symmetries of the
 !! simplex). The corresponding permutation tables are stored in
 !! <tt>nodes\%stab</tt>. The nodes need not to be equally spaced, and
 !! Fekete points could be used. However, for the present
 !! implementation equally spaced nodes are used for simplicity. One
 !! should not, however, rely on this fact, since it should be
 !! possible to switch to Fekete points without changing anything in
 !! the rest of the code in the future. Concerning the order used to
 !! store the nodes in the output array \c nodes, see the comments
 !! about \c t_lagnodes.
 !!
 !! Concerning the implementation details, we note the following.
 !! <ul>
 !!  <li> All the computations are carried out in integer arithmetic
 !!  to allow comparing two nodes for equality, and then the nodes are
 !!  normalized on \f$[0,1]\f$.
 !!  <li> Integer coordinates vary in \f$\{0,\,1,\,\ldots,\,k\}\f$.
 !!  <li> the complete set of nodes for a \f$d\f$-simplex is given by
 !!  all the \f$d\f$-tuples \f$(x_1,\ldots,x_d)\f$ for
 !!  \f{displaymath}{
 !!    x_i\geq 0, \qquad \sum_{i=1}^dx_i \leq k,
 !!  \f}
 !!  while the internal nodes are those for which
 !!  \f{displaymath}{
 !!    x_i\geq 1, \qquad \sum_{i=1}^dx_i \leq k-1,
 !!  \f}
 !!  see also the exponents of the monomial basis computed by \c
 !!  all_sympoly in \c mod_sympoly.
 !!  <li> Introducing an additional coordinate
 !!  \f$x_0=k-\sum_{i=1}^dx_d\f$, a permutation of the \f$d+1\f$
 !!  vertices of \f$\hat{K}\f$ results in the same permutation of
 !!  the coordinates \f$x_0,x_1,\ldots,x_d\f$ (barycentric
 !!  coordinates).
 !!  <li> Starting from the nodes with coordinates such that
 !!  \f$x_i\geq x_{i+1}\f$, all the remaining nodes can be generated
 !!  with permutations. This is useful because the permutation table
 !!  is then obtained as a byproduct.
 !!  <li> Let us consider the generation of the nodes with \f$x_i\geq
 !!  x_{i+1}\f$. Assuming that \f$x_0,\ldots,x_{i-1}\f$ are already
 !!  fixed and letting \f$s_i=\sum_{j=0}^{i-1}x_j\f$, the
 !!  normalization condition implies
 !!  \f{displaymath}{
 !!    x_i = k-s_i-\sum_{j=i+1}^dx_j.
 !!  \f}
 !!  The maximum and minimum values are obtained as follows:
 !!  <ul>
 !!   <li> \f$x_i^{\max}\f$ is obtained for \f$x_j\f$ minimum, i.e.
 !!   \f$x_j=1\f$, yelding
 !!    \f{displaymath}{
 !!      x_i^{\max} = \min(x_{i-1},k-s_i-(d-i)),
 !!    \f}
 !!   <li> \f$x_i^{\min}\f$ is obtained for \f$x_j\f$ maximum, i.e.
 !!   \f$x_j=x_i^{\min}\f$, yelding
 !!    \f{displaymath}{
 !!      x_i^{\min} = \frac{k-s_i}{d+1-i}.
 !!    \f}
 !!  </ul>
 !!  The fact that \f$x_0\geq1\f$, together with the normalization
 !!  condition, also implies the last constraint
 !!  \f$\sum_{i=1}^dx_i\leq k-1\f$.
 !! </ul>
 !<
 pure subroutine lag_nodes(nodes,d,k)
  integer, intent(in) :: d, k
  type(t_lagnodes), allocatable, intent(out) :: nodes(:)

  logical :: node_present
  integer :: id, in, ip, inn, nn, nn_old, ptot
  integer, allocatable :: x(:,:), pcoord(:), coord(:,:,:), &
    perm2coord(:,:), coord2perm(:,:), gen_nodes(:)
  type(t_perm), allocatable :: pi_tab(:)

   allocate( nodes(1:d) )
   do id=1,d
     nodes(id)%d = id

     ! first compute the basic nodes
     call basic_nodes(id,k,x)

     ! now generate additional nodes by permutations
     call perm_table(pi_tab,id+1) ! vertex permutations
     allocate( pcoord( 0:id                          ) , &
                coord( 0:id,size(pi_tab),size(x,2)) , &
           perm2coord(      size(pi_tab),size(x,2)) , &
           coord2perm(      size(pi_tab),size(x,2)) , &
            gen_nodes(                      size(x,2)) )

     do in=1,size(x,2) ! loop over the nodes
       nn = 0 ! counter of the generated nodes
       do ip=1,size(pi_tab) ! loop over the permutations
         ! careful: indexes in x start from 0
         pcoord = x(pi_tab(ip)%i-1,in) ! permuted coords
         ! check whether the node is already present
         node_present = .false.
         check: do inn=1,nn
           if(all(coord(:,inn,in).eq.pcoord)) then
             node_present = .true.
             nn_old = inn
             exit check
           endif
         enddo check
         if(.not.node_present) then
           nn = nn+1
           coord(:,nn,in) = pcoord
           coord2perm(nn,in) = ip
           perm2coord(ip,in) = nn ! how the node was generated
         else
           perm2coord(ip,in) = nn_old ! how the node was generated
         endif
       enddo
       gen_nodes(in) = nn ! number of nodes generated by in
     enddo

     nodes(id)%nn = sum(gen_nodes)
     allocate(nodes(id)%x(nodes(id)%d,nodes(id)%nn))
     allocate(nodes(id)%stab(size(pi_tab),nodes(id)%nn))
     nn = 0
     do in=1,size(x,2)
       do inn=1,gen_nodes(in)
         nn = nn+1
         nodes(id)%x(:,nn) = real(coord(1:id,inn,in),wp)/real(k,wp)
         do ip=1,size(pi_tab)
           ! We know that coord is obtained from x with permutation
           ! coord2perm(inn,in); if we now apply pi_tab(ip) to coord
           ! we get the same result as if applying ptot to x.
           ptot = idx(pi_tab(ip)*pi_tab(coord2perm(inn,in)))
           nodes(id)%stab(ip,nn) = sum(gen_nodes(1:in-1)) + &
                                   perm2coord(ptot,in)
         enddo
       enddo
     enddo

     deallocate( pcoord , coord , perm2coord , coord2perm , gen_nodes )
     deallocate(pi_tab)

   enddo

 end subroutine lag_nodes

!-----------------------------------------------------------------------

 pure subroutine basic_nodes(d,k,x)
  integer, intent(in) :: d, k ! dimension and order
  integer, allocatable, intent(out) :: x(:,:)

  integer :: nxi, m, i, j, si, cjs, cje
  integer, allocatable :: ximin(:), ximax(:), nxi_new(:), temp(:,:)

   ! first coordinate x_0
   allocate(ximin(1),ximax(1))
   ximin(1) = ceiling( real(k,wp)/real(d+1,wp) )
   ximax(1) = k-d
   nxi = ximax(1)-ximin(1)+1

   ! if there are no nodes exit immediately
   if(nxi.lt.1) then
     allocate(x(0:d,0))
     return
   endif
     
   allocate(x(0:0,nxi))
   x(0,:) = (/ (m, m=ximin(1),ximax(1)) /)
   deallocate(ximin,ximax)
   
   ! other coordinates x_i, i=1,...,d
   do i=1,d-1

     ! count the new values
     allocate(ximin(nxi),ximax(nxi),nxi_new(nxi))
     do j=1,nxi
       si = sum(x(0:i-1,j))
       ximin(j) = ceiling( real(k-si,wp)/real(d+1-i,wp) )
       ximax(j) = min( k-si - (d-i) , x(i-1,j) )
       nxi_new(j) = ximax(j)-ximin(j)+1
     enddo

     ! store a temp. copy and reallocate x
     allocate(temp(size(x,1),size(x,2))); temp = x; deallocate(x)
     allocate(x(0:i,sum(nxi_new)))

     ! set the updated values of x
     do j=1,nxi
       cjs = sum(nxi_new(1:j-1)) + 1
       cje = sum(nxi_new(1:j  ))
       do m=cjs,cje
         x(0:i-1,m) = temp(:,j)
       enddo
       x(i,cjs:cje) = (/ (m, m=ximin(j),ximax(j)) /)
     enddo
     nxi = sum(nxi_new)
     deallocate(ximin,ximax,nxi_new,temp)

   enddo

   ! add the last coordinate
   if(d.ge.1) then
     allocate(temp(size(x,1),size(x,2))); temp = x; deallocate(x)
     allocate(x(0:d,size(temp,2))); x(0:d-1,:) = temp
     x(d,:) = k - sum(temp,1)
     deallocate(temp)
   endif

 end subroutine basic_nodes

!-----------------------------------------------------------------------
 
 !> This subroutine can be used to add some new nodes to an existing
 !! set of nodes. The coordinates of the new nodes are passed in \c
 !! new_x.
 !!
 !! If the number of rows of \c new_x is different from
 !! <tt>nodes\%d</tt> the output is an unallocated object,
 !! corresponding to an error.
 !!
 !! The input argument \c stab should be used if the new nodes
 !! represent a closed set with respect to the symmetries of the
 !! elements, and in such a case of \c stab will be corrected to
 !! account for the renumbering which takes place when inserting the
 !! new nodes in the original set. If \c stab is absent, the
 !! corresponding field <tt>nodes\%stab</tt> will contain <tt>-1</tt>
 !! in the columns corresponding to the new nodes.
 !<
 pure subroutine add_nodes(nodes,new_x,stab)
  type(t_lagnodes), intent(inout) :: nodes
  real(wp), intent(in) :: new_x(:,:)
  integer, intent(in), optional :: stab(:,:)

  integer :: stab_temp(size(nodes%stab,1),nodes%nn+size(new_x,2))
  real(wp) ::   x_temp(           nodes%d,nodes%nn+size(new_x,2))

   if(size(new_x,1).ne.nodes%d) then ! invalid
     nodes%d = -1
     nodes%nn = -1
     deallocate(nodes%x,nodes%stab)
     return
   endif

   x_temp(:,1:nodes%nn  ) = nodes%x
   x_temp(:,nodes%nn+1: ) = new_x

   stab_temp(:,1:nodes%nn) = nodes%stab
   if(present(stab)) then
     stab_temp(:,nodes%nn+1: ) = stab + nodes%nn
   else
     stab_temp(:,nodes%nn+1: ) = -1
   endif

   nodes%nn = size(x_temp,2)

   deallocate(nodes%x)
   allocate(nodes%x(size(x_temp,1),size(x_temp,2)))
   nodes%x = x_temp

   deallocate(nodes%stab)
   allocate(nodes%stab(size(stab_temp,1),size(stab_temp,2)))
   nodes%stab = stab_temp

 end subroutine add_nodes
 
!-----------------------------------------------------------------------

 subroutine write_me(me,var_name,fu)
 ! octave output for the master element object
  integer, intent(in) :: fu
  type(t_me), intent(in) :: me
  character(len=*), intent(in) :: var_name

  character(len=*), parameter :: &
    this_sub_name = 'write_me'
 
   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 9' ! number of fields

   ! field 01 : d
   write(fu,'(a)')      '# name: d'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(me%d,'<cell-element>',fu)

   ! field 02 : isv
   write(fu,'(a)')      '# name: isv'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(me%isv,'<cell-element>',fu)

   ! field 03 : children
   write(fu,'(a)')      '# name: children'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(me%children,'<cell-element>',fu)

   ! field 04 : pi_tab
   write(fu,'(a)')      '# name: pi_tab'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(me%pi_tab,'<cell-element>',fu)

   ! field 05 : xv
   write(fu,'(a)')      '# name: xv'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(me%xv,'<cell-element>',fu)

   ! field 06 : n
   write(fu,'(a)')      '# name: n'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(me%n,'<cell-element>',fu)

   ! field 07 : a
   write(fu,'(a)')      '# name: a'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(me%a,'r','<cell-element>',fu)

   ! field 08 : vol
   write(fu,'(a)')      '# name: vol'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(me%vol,'<cell-element>',fu)

   ! field 09 : voldm1
   write(fu,'(a)')      '# name: voldm1'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(me%voldm1,'<cell-element>',fu)

 end subroutine write_me

!-----------------------------------------------------------------------

 subroutine write_simpl_list(simpl_list,var_name,fu)
 ! octave output for t_simpl_list objects
  integer, intent(in) :: fu
  type(t_simpl_list), intent(in) :: simpl_list(:)
  character(len=*), intent(in) :: var_name

  integer :: i
  character(len=*), parameter :: &
    this_sub_name = 'write_simpl_list'
 
   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 3' ! number of fields

   ! field 01 : d
   write(fu,'(a)')      '# name: d'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(simpl_list)
   do i=1,size(simpl_list)
     call write_octave(simpl_list(i)%d,'<cell-element>',fu)
   enddo

   ! field 02 : ns
   write(fu,'(a)')      '# name: ns'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(simpl_list)
   do i=1,size(simpl_list)
     call write_octave(simpl_list(i)%ns,'<cell-element>',fu)
   enddo

   ! field 03 : s
   write(fu,'(a)')      '# name: s'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(simpl_list)
   do i=1,size(simpl_list)
     call write_octave(simpl_list(i)%s,'<cell-element>',fu)
   enddo

 end subroutine write_simpl_list

!-----------------------------------------------------------------------

 subroutine write_lagnodes(nodes,var_name,fu)
 ! octave output for t_lagnodes objects
  integer, intent(in) :: fu
  type(t_lagnodes), intent(in) :: nodes(:)
  character(len=*), intent(in) :: var_name

  integer :: i
  character(len=*), parameter :: &
    this_sub_name = 'write_lagnodes'
 
   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 4' ! number of fields

   ! field 01 : d
   write(fu,'(a)')      '# name: d'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(nodes)
   do i=1,size(nodes)
     call write_octave(nodes(i)%d,'<cell-element>',fu)
   enddo

   ! field 02 : nn
   write(fu,'(a)')      '# name: nn'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(nodes)
   do i=1,size(nodes)
     call write_octave(nodes(i)%nn,'<cell-element>',fu)
   enddo

   ! field 03 : x
   write(fu,'(a)')      '# name: x'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(nodes)
   do i=1,size(nodes)
     call write_octave(nodes(i)%x,'<cell-element>',fu)
   enddo

   ! field 04 : stab
   write(fu,'(a)')      '# name: stab'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(nodes)
   do i=1,size(nodes)
     call write_octave(nodes(i)%stab,'<cell-element>',fu)
   enddo

 end subroutine write_lagnodes

!-----------------------------------------------------------------------

end module mod_master_el

